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Modeling of Turbulent Free Shear Flows 


Dennis A. Yoder, James R. DeBonis, and Nicholas J. Georgiadis 
National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

The modeling of turbulent free shear flows is crucial to the simulation of many aerospace 
applications, yet often receives less attention than the modeling of wall boundary layers. 
Thus, while turbulence model development in general has proceeded very slowly in the past 
twenty years, progress for free shear flows has been even more so. This paper highlights 
some of the fundamental issues in modeling free shear flows for propulsion applications, 
presents a review of past modeling efforts, and identifies areas where further research is 
needed. Among the topics discussed are differences between planar and axisymmetric 
flows, development versus self-similar regions, the effect of compressibility and the evolu- 
tion of compressibility corrections, the effect of temperature on jets, and the significance 
of turbulent Prandtl and Schmidt numbers for reacting shear flows. Large eddy simula- 
tion greatly reduces the amount of empiricism in the physical modeling, but is sensitive 
to a number of numerical issues. This paper includes an overview of the importance of 
numerical scheme, mesh resolution, boundary treatment, sub-grid modeling, and filtering 
in conducting a successful simulation. 
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x c potential core length shifted to align peak centerline turbulence values 

Xi spatial coordinates 

xw potential core length given by Witze relation 
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Subscripts / Superscripts 

D deviatoric tensor, S® = Sij — \Skk$ij 

L laminar 

T turbulent 

sgs sub-grid scale 

1, 2 mixing layer streams 

Operators 

{} trace of the contained tensor expression 

u RANS time averaged or LES spatial filtered velocity 

u RANS density-weighted time average or LES density-weighted spatial filtered velocity 

v! RANS fluctuating velcotiy 

u" RANS density-weighted Favre-fluctuating velocity 


I. Introduction 

W hen asked what are the most challenging turbulent flow problems facing computational fluid dynamics 
(CFD) for aerospace applications, those involved in model development and application would probably 
identify one of many wall bounded issues. These include the difficulty in predicting adverse pressure gradient 
flows, flow separation, and reattachment; streamline curvature; corner flows; shock wave boundary layer 
interaction; turbulence transition; and heat transfer. Each of these contribute to the ability to predict 
aerodynamic drag, engine inlet performance metrics such as pressure recovery, and thermal loading. However, 
turbulent free shear flows such as mixing layers, jets, and wakes also play an important role in aeropropulsion 
applications. These types of flow involve the motion of fluid that is away from solid surfaces. 

Modern subsonic commercial aircraft are typically configured with high-bypass ratio turbofan engines in 
which the bypass fan stream mixes with the high energy exhaust flow of the core engine. This mixing may 
occur in the aft portion of the nozzle or in the plume region. The level of noise generated by these exhaust 
flows increases non-linearly with jet velocity and correlates with the turbulent kinetic energy in the shear 
layer. Many concepts for reducing jet noise have been investigated over the past several years and resulted 
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in significant reductions in jet noise relative to two decades ago. Lobed mixers, chevron nozzles, tabs, fluidic 
injection, and plasma actuators are some of the concepts that have been explored to increase mixing of the 
exhaust streams and modify the turbulence characteristics in the jet. Other concepts have been designed to 
modify the directivity of the noise away from ground observers. This is accomplished by offsetting the core 
and bypass streams or inserting flow deflection devices into the bypass stream in order to divert more of the 
lower-energy fan flow below the core stream. 

With regulations on aircraft noise becoming ever more stringent, the ability to predict and potentially 
modify the flowfield to eliminate these sources has been an active research area. Tools based upon an acoustic 
analogy use the mean flow and turbulence fields from CFD solutions of the Reynolds-averaged Navier-Stokes 
(RANS) equations to estimate noise levels. Alternatively, large-eddy simulation (LES) can be used to directly 
compute the unsteady pressure fluctuations in the jet. In these simulations, prediction of both the mean 
flow and turbulence fields are needed in order to assess the noise of the exhaust nozzle. 

Knowledge of the spreading rate of the exhaust flow is also needed in order to address potential plume 
interaction with the control surfaces of the aircraft and to assess the vulnerability of military aircraft to 
infrared plume signature detection. For air-breathing hypersonic vehicles, the size and weight of the super- 
sonic combustion ramjet (scramjet) are driven by how rapidly the fuel and oxidizer can be mixed to enable 
complete combustion. 

Wake flows are another important type of free shear flow. As aircraft approach for landing, the unsteady 
wake generated by the landing gear and control surfaces of the wing generate considerable noise. As with jets, 
wakes can also have undesirable effects on the downstream control surfaces of the vehicle. In turbomachinery 
applications, the close proximity of the rotor/stator stages can lead to blade row interaction as the wake 
from one blade row impacts the next one downstream. 

In all of these applications, computational fluid dynamics is a critical tool to be used in evaluating future 
concepts. In order to predict the viability and efficiency of these systems one must be able to accurately 
predict turbulent free shear flows. However, simulation of these flows remains a challenge for current CFD 
methods. RANS-based techniques have been used to reproduce trends that occur due to changing geometry 
or flow conditions, but their accuracy is still lacking. This is due, in part, to the fact that most RANS 
models have been developed and calibrated for highly idealized free shear flows and/or the prediction of 
boundary layers. Modifications intended to correct certain model deficiencies have largely been empirical 
in nature, not broadly applicable, and in some cases based upon physical arguments that have since been 
shown to be incorrect. LES and hybrid RANS/LES techniques have demonstrated promise for improving the 
predictive capability for free shear flows. However, there are still a number of uncertainties in such methods 
that prevent their use on a regular basis by the aerodynamics community. 

This paper reviews the current physical modeling capability for two basic free shear flows, namely mixing 
layers and jets, which are similar in nature. These “simple” flows involve more than just a mean shear and 
often display well-organized turbulent structures. The discussion begins with the theory and experimental 
observations of canonical incompressible mixing layers and jets, followed by the changes that are known 
to occur in higher-speed compressible flows more applicable to propulsion applications. Specific modeling 
challenges for both RANS and LES are reviewed, and some suggestions for future research are made. 


I. A. Mixing Layers 


Planar mixing layers form when two parallel streams initially separated by a solid surface come into contact. 
Velocity differences between the two streams create a shear layer, which grows in thickness with downstream 
distance. Far downstream the flow reaches a state where the profiles of mean velocity and turbulence 
statistics achieve self-similarity when plotted using appropriately scaled variables. Generally, the mean 
velocity becomes self-similar before the turbulence quantities, which develop progressively slower in the 
order of streamwise turbulence intensity, transverse turbulence intensity, then turbulent shear stress. 1 The 
shape of the experimentally measured mean velocity profile is often compared with the theoretical error 
function profile derived by Gortler. 2 

u* = - — — = \ [1 + erf (17-770)] (1) 

where the similarity variable rj is defined as 


77 = cr 


y 

X — Xq 


(2) 
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To use this expression, three parameters must be known or determined experimentally: a, the spread rate 
parameter; Xq, the location of the virtual origin of the mixing layer, which is often upstream of the physical 
location where the mixing layer begins; and y o (or r)o), the transverse location of the center of the shear 
layer. The spread rate parameter can be related to the velocity ratio (r) and density ratio (s) of the two 
streams, 3 


(To = , = (1 ~ r) (1 + y / s) 

a s 2(l + r v / s) 


( 3 ) 


and a value of cto = 11 is often used. 

From a practical perspective, there is no way to measure the exact width of the unsteady shear layer. 
Therefore, researchers have devised various measures of the shear layer width which are readily quantifiable 
and can be related to the the spread parameter. One measure of the shear layer width is the visual thickness 
S v i s obtained from images of the flow. The spreading angles of both sides of the mixing layer are estimated 
by drawing straight-line mean tangents to the “edges” of the mixing layer. The visual thickness at any given 
axial location is then defined as the distance between the edge lines. As a more objective measure, researchers 
have used thickness definitions based upon profiles of pitot pressure or mean velocity. The various velocity- 
based thickness definitions include: the 10% Art thickness, 6; the vorticity thickness, (5 U ; the momentum 
thickness, 6\ and the energy thickness, B. Samimy and Elliot 4 discuss some of the benefits and shortcomings 
of these various measures of the mixing layer thickness, particularly with regards to: (1) the number and 
reliability of data points used in describing the thickness, and (2) the ability of that parameter to describe 
variations in the mixing layer thickness due to changes in flow characteristics. 

In the fully-developed region, the mixing layer thickness is found to grow linearly with downstream 
distance. Using equations (1-3), the growth rates can be related to the velocity-density parameter, A s . 


db/dx d5 u /dx dO/dx dB/dx . . 

0.165 - 0.161 - 0.036 “ 0.135 s ^ 

Figure 1 shows that Eq. (4) predicts the growth rate fairly well, except for the last few data sets that involve 
streams with differing densities. Empirical relations for the mixing layer growth rate have also been put 
forth. 3 ’ 5 


dSpa /dx dS V i S / dx d6 u /dx 
0.28 = 0.34 = 0.17 


( 5 ) 


These relations are based upon experimental data rather than Gortler’s velocity profile. 


I.A.l. Experimental Data 

Similarity profiles of mean velocity and turbulence quantities from a number of incompressible mixing layer 
experiments are presented in figure 2. The shape of the experimentally measured mean velocity profiles 
compares very well with Gortler’s error function profile. The turbulence profiles all exhibit very similar 
Gaussian-like curves, which are in line with that obtained analytically for the shear stress when a simple 
eddy viscosity model is used. 6,32 


— ^ = ~P7= — ex P \~ (d ~ ? ?o) 2 l 
(AC/) 2 40ta o L J 

Peak values of the normalized turbulence quantities are found to lie within the approximate ranges: 


0.010 < 

— uV/(A U) 2 

< 0.013 

0.025 < 

«V/(AF) 2 

< 0.034 

0.016 < 

uV/(A U) 2 

< 0.020 

0.020 < 

w'w'/{AU) 2 

< 0.022 

0.032 < 

k/(AU) 2 

< 0.035 


( 6 ) 


( 7 ) 


According to the data of Mehta 25 the normalized turbulence quantities exhibit no discernible trends with 
increasing velocity ratio. However, elevated freestream turbulence levels, such as in the case of Pui and 
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Figure 1: Incompressible mixing layer growth rates from reference 6. Solid symbols and labels indicate 
density ratios other than 1.0. 


Gartshore, 20 can lead to higher peak turbulence values within the mixing layer and an increase in the mixing 
layer growth rate. 23,32 References 35 and 24 investigated the effect of the boundary layer state upstream of 
the mixing layer and found that when laminar boundary layers are present the growth rate and Reynolds 
stresses are significantly higher in the development region than when the boundary layers are turbulent. 
Far downstream the turbulence values for both cases asymptote to approximately the same constant values. 
Champagne 1 ' conjectures that if the initial boundary layer is turbulent, then the orderly vortex-pairing 
process (described in the next section) might occur randomly, in some intermittent fashion, or not at all, 
and therefore lead to subsequent differences within the development region. 

I. A. 2. Structure 

In incompressible turbulent mixing layers, large-scale spanwise-oriented vortical structures have been ob- 
served. 5 The formation of these large-scale structures is triggered by the fundamental Kelvin-Helmholtz 
instability of the flow, and it is the exponential growth of these instabilities which causes the shear layer to 
roll-up into large spanwise vortices. The initial spacing between vortices is tied to the frequency of the most 
dominant instabilities, as evidenced when forcing is applied to the mixing layer. 

These spanwise rollers convect downstream at roughly the average speed of the two streams. As they 
do, they grow by entraining fluid from outside the shear layer and they combine through a series of pairing 
and tearing processes. 36 In the pairing process, two adjacent eddies draw together and begin rotating about 
each other, redistributing their vorticity until they eventually merge into a single larger eddy. This pairing 
process may be triggered by local instabilities which push one of the eddies away from the centerline, thereby 
giving it a slightly faster or slower convective velocity depending on which stream it is pushed towards. In the 
tearing process, a portion of an eddy splits apart from the rest and is consumed by another neighboring eddy. 
This amalgamation of eddies through the pairing and tearing processes results in fewer, larger structures. 
Therefore, the spacing between eddies increases with distance downstream. 

As the flow progresses downstream, secondary spanwise instabilities are also triggered which lead to 
the development of additional streamwise vortices. 37 These counter-rotating streamwise vortices or “ribs” 
form in the braid region between rollers and wrap around successive spanwise vortices. 38,39 The spanwise 
spacing between these streamwise vortices is observed to increase further downstream and may be related 
to the increasing streamwise spacing of the larger rolling structures. 40 Jimenez, Cogollos, and Bernal 39 
estimate that the strength of the streamwise ribs can be a significant fraction of the local circulation in the 
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k/(U,-U 2 ) z -uV/(U,-U 2 ) 2 U' = (U-U 2 ) / (U,-U 2 ) 
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Eq. (6) 
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(c) Turbulent Kinetic Energy 

Figure 2: Similarity profiles of incompressible mixing layers from reference 6. 
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(f) Spanwise Normal Stress 

Figure 2: Similarity profiles of incompressible mixing layers from reference 6 (continued). 
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spanwise eddies. The interaction of these structures causes a distortion of the large-scale spanwise rollers, 
three-dimensionality, and vortex stretching. These factors work to feed higher-order instabilities which 
lead to the development of smaller scale turbulence. The energy of the large turbulent eddies cascades to 
progressively smaller eddies until it reaches the smallest Kolmogorov scales where it is then dissipated into 
heat by molecular viscosity. This action leads to the eventual decay of all of the turbulent structures. 


I. A. 3. Compressibility 


One of the key observations of compressible planar mixing layers is that they tend to spread more slowly than 
incompressible mixing layers at equivalent velocity and density ratios. Since many of the early experiments 
investigated single-stream mixing layers, where one stream was at rest and increasing the Mach number of 
the other was accompanied by decreasing temperature and therefore increasing density, this reduction in 
growth rate was generally thought to be caused by density differences between the two streams. However, 
subsequent dual-stream experimental results 5 revealed that the effect of extreme density differences alone 
could not explain this phenomenon and that the reduction in growth rate of compressible mixing layers is a 
true compressibility effect. 

The convective Mach number is one of the most commonly used measures for quantifying the level of 
compressibility of mixing layer flows. For two pressure-matched streams with the same ratio of specific heats, 
the convective velocity and Mach number are computed as 


U c = 


a 2 u\ + a\u 2 


M r = 


Cl 1 + Ct2 

Hi - U 2 

ai + a 2 


( 8 ) 

(9) 


and represent the rate at which the large-scale structures propagate downstream. 

Using the convective Mach number as a measure of compressibility, the mixing layer growth rates from 
a variety of experimental data sets are summarized as shown in figure 3. These growth rates have been 
normalized by the correlation for incompressible mixing layer growth rate at the same velocity and density 
ratios. The figure indicates a clear pattern of reduced growth rate with increasing compressibility, at least 
until M c > 1.0 where additional compressibility seems to have only a minimal effect. The scatter in the 
data has been attributed to a number of factors, including: the use of different experimental facilities and 
measurement techniques, the use of data which has not fully achieved self-similarity, differences in the mixing 
layer thickness definition used, uncertainty in the estimate of the incompressible mixing layer thickness, and 
the possibility that the convective Mach number may not fully capture the effects of compressibility or may 
only capture the first-order effects. 38 In particular, the data of figure 3 indicate a disparity between mixing 
layer growth rate measurements based on pitot thickness or visual thickness (shown as solid symbols) and 
the various other growth rate estimates which are based on velocity measurements. 
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Figure 3: Normalized compressible mixing layer growth rate data from reference 6. 


Figure 4 compares available Reynolds stress data for compressible mixing layers. There is a significant 
reduction in the shear stress and transverse normal stress with increasing convective Mach number. These 
changes in turbulence statistics occur at convective Mach numbers as low as 0.2, well before any reduction 
in the corresponding mixing layer growth rate data is observed (M c « 0.5). Trends in the data for the 
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streamwise turbulence intensity are not as straightforward. The measurements of Goebel and Dutton 1 and 
Gruber, et al. 48 indicate that the streamwise turbulence remains relatively unaffected by compressibility. 
Others have noted a slight reduction, though not as severe as with the other turbulence quantities. This 
disparity in the measured change in turbulence anisotropy is quite evident by plotting the ratio u'u r /v'v' 
versus convective Mach number, as done in figure 4d. 
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Figure 4: Turbulent stress values for compressible mixing layers from reference 6. 


Along with a reduction in mixing layer growth rate and turbulence statistics, structural changes have 
also been observed in compressible mixing layers. 58,59 At low convective Mach numbers, large-scale vortical 
structures are sometimes observed to extend spanwise across the width of the test section much like they 
do in incompressible flows. At moderate convective Mach numbers (M c « 0.5) the Brown-Roslrko style 
large-scale organized structures are less frequently encountered and often appear smaller in size. Instead 
the high- and low-speed interfaces of the shear layer appear less smooth with the freestream fluid sharply 
intruding into the middle of the shear layer. Plan views indicate that the mixing layer becomes more three- 
dimensional with less observed spatial regularity. At higher convective Mach numbers, the structures are 
highly three-dimensional and the flows tend to exhibit thin streamwise sheets of vorticity which appear at 
various transverse locations. These sheets do not extend as far across the mixing layer as the incompressible 
rolling structures and exhibit little transverse communication among themselves. It has been suggested that 
this reduced communication 60,61 across the shear layer results in the reduction of the transverse turbulence 
intensity and turbulent shear stress, which affects the mixing layer growth rate. 
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I.B. Jets 


In jet flows, one fluid stream mixes with a surrounding medium that is either at rest or in motion. Viscous 
forces along the shared interface between the inner and outer streams result in the formation of a shear layer 
that transfers momentum between streams. The largely inviscid flow within the inner stream is referred to 
as the potential core. The shear layer grows with distance downstream until it completely consumes the 
potential core. For jets issuing from the round nozzles found on many aircraft engines or rocket exhaust 
systems, the shear layer is circular in shape. Flow in the initial part of the jet may be loosely characterized 
as a planar-like shear layer, because the width of the shear layer is small compared to the diameter of the jet. 
Gutmark, Schadow, and Yu 62 discuss some of the differences between planar mixing layers and jet flows. In 
particular they note that round jets support additional axisymmetric (ring-like) and helical modes and that 
the number of vortex interactions in jets is limited by the introduction of additional length scales, such as 
the length of the potential core. Linearized stability analysis has demonstrated that both the axisymmetric 
and helical modes of jets become more stable as the Mach number increases. 63 

Jets from non-circular nozzles support a different set of instability modes which are parametrized by 
factors such as eccentricity, aspect ratio, circumferential distribution of momentum thickness, and local 
radius of curvature. 64 In these flows, shear layer growth and the formation of small-scale turbulence is 
driven more by the dynamics of vortex self-induction, also known as (^-dynamics, than the vortex pairing 
process. 65 Self-induction 66 refers to the phenomenon whereby the advection velocity of a curved vortex 
filament is moves faster in regions of greater curvature. For elliptical or rectangular jets, that portion of the 
azimuthal vortex ring located near either end of the major axis will accelerate out-of-plane. The distortion 
of the vortex ring causes an induced velocity on the minor-axis sides which pushes the flow outward. After 
several contortions, the vortex ring becomes planar once again but appears to be rotated by 90 degrees 
relative to the initial orientation. This axis switching can repeat several times before the vortex breaks 
downs and decays. 

For nozzles that have corners or flow control devices such as tabs, the induced motion of streamwise 
vortex pairs in the flow, known as w x -dynamics, is important. 67 In rectangular nozzles, differences in the 
cross-stream pressure gradients inside the nozzle result in the formation of secondary flows consisting of 
streamwise-oriented vortex pairs. The induced motion of these vortex pairs will cause them to move from 
the corners toward the center of the jet and displace flow towards the sides. This dynamic can likewise lead 
to an axis switchover, but will not support further switchings beyond the first. Nozzles that have both flat 
sides and corners benefit from the large-scale mixing that occurs at the flat sides and the small-scale mixing 
due to the corner flow. 

The reduced shear layer growth rate in compressible turbulent jets results in longer potential core lengths. 
In practice, jets are found to exhibit compressibility effects at jet Mach numbers higher than about 0.5. For 
unheated incompressible jets the potential core length has been found to be approximately 5 nozzle exhaust 
diameters. For an unheated Mach 2 jet, the potential core length is approximately 10 nozzle diameters long, 
twice the value for an incompressible jet. 

It has been experimentally observed that for otherwise similar jets, e.g. from the same nozzle at same 
jet exit Mach number, a heated jet will exhibit faster mixing than the corresponding unheated jet. Lau 68 
and Lepicovsky 69 concluded that the boundary layer state at the end of nozzle my be largely responsible 
for the differences in heated jet development. However, other studies such as that of Monkewitz and Sohn'° 
and Monkewitz, et al.' 1 have shown that the instabilities in the initial part of a heated jet give rise to a 
faster rate of turbulent structure growth and resultant faster jet mixing rate in comparison to unheated jets. 
They indicate that for jets with temperature ratios greater than 1.4 (or corresponding density ratios less 
than approximately 0.7), a local absolute instability region may develop in the jet potential core which is 
primarily responsible for this enhanced mixing behavior. 

Witze 72 examined a number of jet experiments to develop expressions for the centerline velocity decay 
and jet potential core length as a function of jet Mach number and ratio of jet to freestream densities. 


U 

Ujet 


= 1 — exp 

Xw 

D 


1 

L0.70 (1 — ar/arw) J 

0.70 


W 2 Poo Ip jet 
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Table 1: Test conditions used in reference 73. 


Set Point 

M a 

Tjet/Too 

NPR 

Mj et 

3 

0.500 

0.950 

1.197 

0.513 

7 

0.900 

0.835 

1.861 

0.985 

23 

0.500 

1.764 

1.102 

0.376 

27 

0.900 

1.764 

1.357 

0.678 

29 

1.330 

1.764 

1.888 

1.001 

46 

0.900 

2.700 

1.219 

0.548 

49 

1.485 

2.700 

1.678 

0.904 


For subsonic jets, 

« = 0.08 (1 - 0.16 M jet ) ( Poo/Pjet )~ °' 22 (12) 

For properly-expanded supersonic jets, 


/ 9 \— 0.15 

k = 0.063 (. Mf et - l) (13) 

This expression is valid from the nozzle exit to the sonic point. Downstream of that, the subsonic expression 
can be used provided an appropriate shift in x is made. 

Bridges and Wernet 73 provide a concise yet thorough review of available data for subsonic round jets. 
They show that normalizing the axial coordinate by Witze’s potential core length leads to a satisfactory 
collapse of low speed axial turbulence data along the jet centerline. They also present their own particle 
image velocimetry (PIV) data for various nozzles operated over the range of conditions listed in table 1. 
After careful consideration of the measurement error, they arrive at a consensus data set that represents the 
best estimate of the flow statistics and uncertainty. 

Figure 5 compares the centerline profiles of mean and turbulence axial velocity for all of the consensus 
data sets. While normalizing by Witze’s correlation for potential core length provides good collapse of most 
of the data, Bridges and Wernet prefer to use a potential core length x c that is deliberately shifted to align 
the peak turbulence values along the centerline. Doing so results in centerline profiles for the streamwise 
turbulence that are similar in shape, but increase in peak amplitude with increasing jet temperature. 



Figure 5: Axisymmetric jet centerline profiles of mean (left) and turbulent (right) axial velocity from reference 
73. 
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Figure 6 shows the sensitivity of both potential core metrics to temperature and Mach number. The 
elevated turbulence levels in heated jets results in shorter potential core lengths. As the jet Mach number 
increases, so do compressibility effects which reduce the jet spread rate and increase the length of the 
potential core. It is also evident that as the jet Mach number increases, so does the discrepancy between the 
two potential core metrics. The core length based on peak centerline turbulence value tends to shift farther 
downstream. 



Figure 6: Axisymmetric jet potential core lengths as a function of Mach number for jets with different 
temperature ratios from reference 73. Squares, unheated; circles, Tj et /T 00 = 1.7; triangles, Tj et /T oc = 2.7. 
Open symbols are x\y/Dj et , closed symbols are x c /Dj et . 


The peak axial turbulent velocity on the centerline and in the entire flowfield is plotted against potential 
core length in figure 7. The peak turbulence values do not occur along the centerline, but rather in the shear 
layer just downstream of the nozzle lip. However the peak centerline value appears to be proportional to the 
overall peak value. Higher turbulence values result in a more rapid growth of the jet shear layer and shorter 
potential core length. 



Figure 7: Axisymmetric jet peak axial turbulent velocity on the centerline and throughout the flowfield as 
a function of potential core length from reference 73. 

Radial spreading of the jet is another important parameter. Mean and turbulent axial velocity contours 
are shown in figure 8. Even though the axial coordinate is scaled by the potential core length, the contours 
show some variation with static temperature (or density). The heated jets do not appear to extend as far 
radially and the centerline velocity decays more quickly. Higher turbulence levels are also observed farther 
downstream and closer to the centerline. Bridges and Wernet show that better collapse of the contours 
can be made if the radial coordinate is scaled by the jet half-width (ro.5) rather than nozzle diameter. The 
half-width, which is a measure of the center of the jet shear layer, is shown to be insensitive to Mach number. 
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However, the center of the shear layer for hot jets is shifted closer to the axis. The shear layer growth rate, 
on the other hand, is sensitive to Mach number but not temperature. For hot jets then, the shear layer is 
shifted closer to the axis but grows toward the centerline at the same rate as cold jets. The result is a shorter 
potential core. 



II. RANS Modeling Issues 


The Reynolds-averaged Navier-Stokes (RANS) equations are routinely used to simulate aeropropulsion 
flows. With this method, a time averaging procedure is used to represent the net effect of the turbulent 
fluctuations without having to simulate the entire unsteady motion. In this respect the method is very 
efficient. However, the averaging process introduces unknown correlation terms that must be modeled. 

The Reynolds averaging procedure decomposes each instantaneous flow variable into a mean and a 
fluctuation about that mean. For example, 

u = u + v! (14) 

where the Reynolds-averaged mean value is defined as 


u = lim 
At — ► » 


1 

At 


to+At 


dt 


(15) 


and At represents a long period of time relative to the fluctuations. For compressible flows, a density-weighted 
procedure leads to a more compact form of the flow equations. Here, 


u = u T u" 


(16) 


where the Favre-averaged mean value is defined as u = puf p. The resultant averaged form of the continuity, 
momentum, and energy equations can be written as: 


dp dpuj 

dt Ox* 


= 0 


dpu. 


dt + d^- t pUiUj + p6ij ~ ( fij + r 5)l : 0 


d 

dt 


p[e + -UiUi + k 


dxj l 


p Uj 


h + m + kj — Ui (jij + r^) + (qj + gj) 


(17) 

(18) 

= 0 (19) 
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where the important correlation terms to be modeled are the turbulent heat flux, the turbulent stress tensor, 
and the turbulent kinetic energy defined as: 


Qj = pu"h" 


T // // 

Tii = ~PU"Uj 


pk = -pu'(u'l 


(20) 

(21) 

(22) 


II. A. Model Calibration 


Turbulence models are approximations used to mimic the mean behavior of the actual physical processes, and 
as such usually involve closure coefficients that must be calibrated through comparison of results with known 
or expected behavior. This is usually done with a series of fundamental flows that isolate specific terms. 
Each turbulence model may have been calibrated using a different procedure, a different set of experimental 
data, or a different metric for “optimizing” the results. Unfortunately, the details of how a model has been 
calibrated are not always fully disclosed. 

To illustrate how models are calibrated, consider the form of a typical k — e eddy viscosity model which 
has been routinely applied to free shear flows. 

Tij = -pu”u” = 2 p, T S? ~\pk Sij (23) 

Pt = Cfjf^pk 2 /e (24) 


d(pk ) d(pujk) 

dt dxj 


d 

dxj 



dk 

dxj 


+ V — pe + Lf~ 


9(pe) d(puje) 
dt dxj 


d 

dxj 



de ' 
dxj _ 


+ Celfl l V 



+ L e 


(25) 

(26) 


V = 


T duj 

Til d^ 


(27) 


In these equations, / M , fi, and /2 represent damping functions and L k and L e represent near-wall terms, the 
details of which are not important for this discussion. The eddy viscosity coefficient C M is often calibrated 
under the assumption of thin shear flows where the production of turbulent kinetic energy is nearly in 
balance with the turbulent dissipation rate ( V « pe) together with the experimental data of Townsend' 4 
which indicates that the ratio of shear stress to turbulent kinetic energy for a wide range of flows has nearly 
the same value of 0.3. The C e 2 coefficient multiplying the dissipation term in the e-transport equation is 
determined by analyzing the decay of isotropic homogeneous turbulence. The C e \ coefficient multiplying the 
production term in the e-transport equation is obtained by applying the model to homogeneous turbulent 
shear flow. The diffusion coefficient oy that also appears in the e-transport equation is calibrated by analyzing 
the log-layer region of the turbulent boundary layer. The diffusion coefficient the fc-transport equation is 
typically set to 1.0. 

An e-based Reynolds stress model requires the calibration of a number of coefficients in the model for 
the pressure-strain correlation tensor. The following quasi-linear form is commonly used. 


n 


D 

ij 

Pt 


~ ( C? + C\ -) bij + C 2 -SP + C 3 - 


ptJ 10 ' e 


bik^kj H - $ikbkj ^buinSmnSij 


Iz r 

(Hi h, n ttk j Rikbkj 


(28) 


The “slow” pressure-strain coefficient C® is calibrated from the decay of anisotropic homogeneous turbu- 
lence. The remaining pressure-strain coefficients Cl, C 3 , and C4 are optimized to provide agreement with 
homogeneous shear flow and rotating shear flow, while C2 is determined from rapid distortion theory for 
homogeneously strained turbulence that is initially isotropic. 
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One shortcoming of any such calibration procedure is that it is assumed that these coefficients will retain 
the same values for all types of flows. Bradshaw 75 notes that the coefficients used in the modeling terms will 
only assume constant values in the unlikely event that the modeled term exactly correlates with the physical 
term being modeled. More plausibly these coefficients will be different for different flows, and may vary from 
point to point within the same flow. He suggests that these coefficients may be improved by making them 
empirical functions of the local flow parameters such as the ratio of production to dissipation or a normalized 
measure of the flow invariants. 

II. B. Differences Between Planar and Axisymmetric Flows 

Free shear flows are composed of a variety of three-dimensional turbulent structures. However, the organi- 
zation of those structures is related to the dominant instability modes, which differ between mixing layers, 
jets, and wakes. These differences can be seen in experimental flow visualization and in solutions from large- 
eddy simulation. This change in the way that the turbulence is organized presents a particular challenge 
for RANS-basecl methods which model statistical averages rather than structure. It is among the reasons 
why most turbulence models cannot accurately predict both mixing layers and jets using the same set of 
coefficients. 

Table 2 compares RANS model predictions for spread rates in the self-similar regions of different types 
of free shear flows. From these results one can see that every turbulence model predicts each type of flow 
with varying degrees of success. For example, the Spalart-Allmaras ' 6 model does a poor job of predicting 
jet flows, but excels in the wake flows for which it was specifically calibrated. In the case of the widely used 
k — e model of Launder and Sharma, ' ' the spread rate is in agreement with the experimental data for the 
mixing layer, plane jet and radial jet, but is too fast for the round jet. 


Table 2: Free shear flow spread rates predicted by RANS models. 



Mixing Layer 

Plane Jet 

Round Jet 

Radial Jet 

Wake 


Experiment 7, 78 81 

0.103-0.120 

0.100-0.110 

0.086-0.096 

0.096-0.110 

0.320-0.400 

k — Ul 

Wilcox 82 * 

0.096 

0.108 

0.094 

0.099 

0.326 

k — Ul 

Wilcox 83 * 

0.105 

0.101 

0.088 

0.099 

0.339 

k — Ul 

Wilcox 84 * 

0.141 

0.135 

0.369 

0.317 

0.496 

SST 

Menter 85t 

0.100 

0.112 

0.127 

- 

0.257 

k — e 

Launder & Sharma 7 '* 

0.098 

0.109 

0.120 

0.094 

0.256 


with Pope 86 * 

0.098 

0.109 

0.086 

0.040 

0.256 

k- c 

Robinson, et al. 8 ' 

0.112 

0.115 

0.091 

0.097 

0.313 

k — t 

Speziale, et al. 88 

0.082 

0.089 

0.102 

0.073 

0.221 

SA 

Spalart & Allmaras 76 * 

0.109 

0.157 

0.248 

0.166 

0.341 


* Results reported by Wilcox . 82 
t Results reported by Bardina, et al . 89 


Pope 86 developed an additional parameter to improve k — e model round jet predictions based on the 
concept of vortex stretching. He theorized that the stretching of turbulent vortex tubes by the mean flow 
has a significant influence on the rate of turbulent dissipation. As these vortices are stretched, there is a 
transfer of energy from the large eddies to the smaller ones. Since the small eddies control the dissipation, the 
dissipation rate e must increase. To achieve this effect, Pope proposed a model correction to the e-transport 
equation that reduces the effective dissipation of e. 

C e2 f 2 p \ - ( c e2 f 2 - C e3 f 3Xp ) (29) 

X P is a nondimensional measure of the vortex stretching. 

X P =(^-J RijRjkSki (30) 

Similar terms have been incorporated into other turbulence models. The Wilcox 82,83 k — u> and Robinson, et 
al. 87 fc — C models include a term of this type into their baseline models. As shown in table 2, these models 
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do reasonably well for each type of free shear flow. Depending upon the behavior of the baseline turbulence 
model, this correction can adversely affect the prediction of other three-dimensional flows where \p 7 ^ 0- 
Such is the case with the k — e model, where the spread rate for the radial jet is significantly reduced. 

For subsonic round jets, baseline k — e models underpredict the growth of the initial shear layer, resulting 
in potential core lengths that are too long. Near the end of the potential core and farther downstream, the 
shear layer growth rate is too high relative to data. Activation of the vortex stretching correction with these 
models will reduce the level of turbulence and slow the shear layer growth rate throughout the flowfield. 
While this action may bring the self-similar mixing rate in line with experimental data, it will have the 
adverse effect of yielding potential core lengths that are even longer than with the baseline model. 

Birch 90 ’ 91 has questioned the physical basis for this correction, noting that experimental results have 
shown lateral divergence actually increases the turbulent shear stress. Regardless, this method is still 
sometimes used because it tends to shift the self-similar RANS predictions in the right direction for round 
jets. 

In algebraic Reynolds stress modeling, the turbulence anisotropy tensor can be related to ten basis tensors 
and six invariants formed from combinations of the mean strain and rotation rate tensors. 92 It is interesting 
to note that the expression RijRjkSki used in Pope’s vortex stretching parameter in Eq. (30) matches one 
of these six invariants: {S' 2 }, {-R 2 }, {S 3 }, {SR 2 }, {S 2 R 2 }, {S 2 R 2 SR}. This lends some support to the 
previously described notion that perhaps the closure coefficients should not be constants, but functions of 
the flow invariants. 

Considering the nature of turbulent jet flows, it is surprising that RANS models do as well as they do. 
The initial plane (or axisymmetric) jet region is similar to planar (or annular) mixing layers. Near the end 
of the potential core, these shear layers interact with each other and merge together. Farther downstream 
the plume resembles a single shear layer. In each region, the composition of turbulent structures differ, and 
one might expect each region to require its own set of coefficients. Defining where and how those coefficients 
should transition would most likely limit the applicability of any such model to other types of flow. 

Over a decade ago, Birch 91 described the dichotomy between Reynolds-averaged methods and structure- 
based methods such as large-eddy simulation. From the coherent structure viewpoint, shear layer growth 
is likely sensitive to the detailed structure of the large eddies. Even though a self-similar region might be 
reached, the turbulent mixing could continue to depend upon the initial turbulent structure and will never 
truly be the same. From the Reynolds-averaged perspective, initial conditions might persist for some distance 
downstream, but an asymptotic state is eventually reached that is unique for each flow geometry (mixing 
layer, jet, etc.). 

Of these two perspectives, it seems more likely that structure-based methods should provide more accurate 
simulations. Large-eddy simulation methods resolve the most important turbulent scales through the use of 
the unsteady Navier-Stokes equations. RANS models are instead based upon a time (or ensemble) average 
of those equations, which as a simplified subset of the full equations will surely contain less information. If 
structure does play such a key role in these types of turbulent flows, then it will be difficult for RANS models 
to accurately replicate such physics even in an averaged sense. 

II. C. The Developing Region 

The types of homogeneous and isotropic benchmark cases used to calibrate the coefficients in RANS-based 
turbulence models are rather idealized and might, at best, resemble that observed in the self-similar region 
that exists far downstream. The developing region is frequently ignored when assessing the performance 
of RANS models, as the spreading rate in the self-similar region is usually emphasized (i.e., table 2). For 
propulsion applications, however, the primary region of interest is often in the near field. As an example, 
the mixing that occurs in scramjets and internal mixer nozzles is almost always done in the minimum length 
possible. Also, most of the noise in a subsonic jet is generated within the approximate length of the potential 
core. 

In the development region, the turbulence is in a non-equilibrium state and could even be transitioning 
from laminar to turbulent. The mean and turbulence profiles vary spatially but, unlike the downstream 
region, are not self-similar. The development region can also be influenced by the upstream flow, including 
inflow profiles, boundary layer state (laminar or turbulent), and freestream turbulence levels. It is unclear 
whether RANS methods can be made sensitive to all of the factors influencing the development region, or 
whether one should use LES instead. 
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II. D. Compressibility 


In RANS simulations, the most common means for including compressibility effects on turbulence is through 
the inclusion of explicit dilatation terms that appear in the transport equations for the Reynolds normal 
stresses (and subsequently the turbulent kinetic energy). These terms, which include the dilatation dis- 
sipation 93 (ed) and the pressure dilatation 94 (n™), have an additional dissipative effect on the growth of 
turbulence. As shown in figure 9 from reference 6, these terms provide varying degrees of growth rate re- 
duction in compressible mixing layers. Both of the dilatation models shown are functions of the turbulent 
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Figure 9: RANS predictions of compressible mixing layer growth rate from reference 6. 


Mach number and display sensitivity to compressibility even at low M c . Therefore, the resulting change in 
growth rate is more gradual than indicated by the data. In addition, these methods do not significantly alter 
the Reynolds stress anisotropy, because the normal stresses are all reduced in an isotropic manner. 

Dembowski and Georgiadis 95 investigated the ability of typical RANS models to predict a Mach 2 round 
jet at different operating temperatures. 96 The results, in figure 10, show that the baseline SST model mixes 
much faster than indicated by the experimental data. This is consistent with the observation that standard 
models predict shear layers that grow too quickly because the models are not sensitive to compressibility 
effects. Use of a dilatation dissipation correction slows the growth of turbulence, but results in a potential 
core length that is too long. Farther downstream the decay rate is still slightly over-predicted for the 
unheated case, but this trend reverses for the higher temperatures. 

Results from direct numerical simulation (DNS) of homogeneous shear flows, 97 temporal mixing layers, 55 
and annular mixing layers 98 indicate that the contributions of dilatation dissipation and pressure dilatation 
are both insignificant. Instead, reduced pressure fluctuations (i.e., reduced communication) across the com- 
pressible shear layer alter the pressure-strain correlation tensor that appears in the Reynolds stress transport 
equation. This tensor redistributes the turbulence energy between the streamwise, transverse, and spanwise 
terms. While each of the Reynolds stresses are reduced with increasing compressibility, the streamwise 
normal stress is less affected and leads to an increase in the turbulence anisotropy. However, for temporal 
mixing layers, the DNS of Pantano and Sarkar 34 suggests that the degree to which turbulence anisotropy 
changes can depend upon the extent and duration of the simulation. 

Since compressibility effects appear to be linked to changes in the pressure-strain correlation tensor, 
accurate physical modeling of compressible mixing layers requires (at a minimum) the use of algebraic or 
differential Reynolds stress models. Formulation of a truly compressible model for the pressure-strain corre- 
lation is a formidable challenge. The more common approach has been to use existing pressure-strain models 
and apply damping functions 34, 99-101 to modify some or all of the pressure-strain model coefficients. These 
damping functions depend upon some measure of compressibility, such as the convective (M c ), turbulent 
(Af t = y/2k/a), or gradient (M g = \S\l/a) Mach number. Of these, the gradient Mach number appears to 
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(a) Centerline axial velocity (b) Centerline stagnation temperature 

Figure 10: Effect of temperature on jet predictions from reference 95. 


be the most appropriate choice since the convective Mach number is a global parameter and the turbulent 
Mach number can be large in both high-speed mixing layers and boundary layers. 

Gomez and Girimaji 102 suggest that modeling of the pressure strain correlation tensor can be split into 
three distinct regimes. At low gradient Mach number, a standard incompressible pressure-strain correlation 
model may be used. At intermediate M g , the pressure and inertial terms are of the same magnitude such 
that the pressure-strain tensor approximately balances production. At high M g , compressibility severely 
limits pressure communication and the contribution of the pressure-strain tensor becomes negligible. The 
transition between regimes is accomplished by a blending function. 

Once again, it is interesting to note that the parameter used to affect a change in the turbulence model 
predictions can be related to the flow invariants. In this case, the gradient Mach number can be related to 
{S 2 }. Thus, procedures that sensitize the model coefficients to M g are, in effect, sensitizing them to one of 
the local flow invariants. 

II. E. Thermal and Mass Effects 

Propulsion flows generally involve gaseous species that are the products of combustion at elevated pressure 
and temperature. As the exhaust expands out the nozzle, it undergoes a significant temperature change. 
Density and thermal gradients in the jet add another level of complication to modeling. 

The review by Georgiadis and DeBonis 103 includes solutions of various RANS models for unheated and 
heated incompressible round jets. Some of the models were specifically tuned for jet flows and/or temperature 
effects. For an unheated jet, baseline two-equation turbulence models generally predict potential core lengths 
that are too long and downstream mixing rates that are too high. While the peak value of turbulent kinetic 
energy along the centerline is well represented, it does not diffuse to the centerline as quickly as the data 
indicates. For the heated jet, the baseline models again predict a potential core that is too long. They also 
predict a centerline turbulence profile that is similar to the unheated case, albeit shifted slightly upstream, 
and fail to predict the higher peak turbulence level. Through a modification of the turbulent viscosity, 
the temperature corrected model mixes more quickly than the standard model, providing better centerline 
velocity decay. However, the turbulent kinetic energy field dissipates much too quickly. 

In RANS simulations, Reynolds’ analogy between momentum and energy is typically used with an as- 
sumed constant turbulent Prandtl number. A value of 0.89 is common for boundary layers, while a value as 
low as 0.5 is more appropriate in free shear flows. Dembowski and Georgiadis 95 investigated the effect of 
turbulent Prandtl number on Mach 2 hot round jet simulations. As shown in figure 11, higher values of Pr T 
have little effect on the velocity field, but slightly reduce the thermal decay. 

Turbulent Prandtl and Schmidt number variations also affect the mixing and combustion that occurs 
in reacting mixing layers. The simulation results 104 for a supersonic reacting mixing layer 105 are shown in 
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Figure 11: Effect of turbulent Prandtl number on hot jet predictions from reference 95. 


figure 12. For this case, the ignition point moves slightly upstream with increasing Pr T . The results show a 
much greater sensitivity to Sc T , which influences the species mass diffusivity much in the same way that the 
turbulent Prandtl number influences the thermal diffusivity. As the turbulent Schmidt number is reduced, 
the species diffuse more quickly creating a more favorable mix of reactants for combustion to occur sooner. 

Unfortunately, the optimal turbulent Prandtl and Schmidt numbers for one case are not optimal for 
others. In general the values will vary throughout the flowfield. Some work has been done in the area of 
variable turbulent Prandtl and Schmidt number modeling. 106-109 Most of these models relate the thermal 
(or mass) diffusivity to the scalar variance of temperature (or species mass fractions), then solve modeled 
transport equations for that variance and its dissipation rate. For combusting flows, development of these 
models is complicated by the consideration of turbulent-chemistry interaction effects, which are sometimes 
included through probability density function (PDF) formulations. A major hindrance to the development 
of variable Pr T and Sc T models is a lack of experimental data to calibrate and validate the thermal and 
mass diffusivities, and underlying scalar transport equations. While these methods are promising, they are 
not yet mature enough for routine use. 

III. LES/DNS Modeling Issues 

Large-eddy simulation and direct numerical simulation hold great promise for improving the predictive 
capability of jets and mixing layers. These scale-resolving methods directly compute the key turbulent 
structures. DNS directly computes the full turbulent spectrum all the way down to the Kolmogorov scales. 
LES uses a spatial filtering procedure to separate the unsteadiness into large, resolved scales and small, 
unresolved scales. LES directly computes the resolved scales and approximates the effect of the unresolved 
scales using a sub-grid model. Because of the very fine grids required to resolve all scales, DNS for most flows 
at realistic Reynolds numbers is impractical using current computer technology. LES significantly relaxes 
the computational requirements by only resolving the large scales and has been used extensively for jets. 
The focus of the discussion here will be on LES techniques. 

In LES, a spatial filter G with a filter width A is used to separate the resolved and unresolved scales. 
For example, 


Ui = 


/ OO 

G(z-0«i(£K (31) 

-OO 

The overbar represents the resolved, filtered, or large-scale portion of the variable. For compressible flows a 
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Figure 12: Effect of turbulent Prandtl number (a,c,e) and Schmidt number 
reference 104. 


(b,d,f) on combustion, from 
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density-weighted filter is often used instead. 


u = pu/p 


(32) 


Applying the filter to the Navier-Stokes equations yields the following forms of the continuity, momentum, 
and energy equations 


d_ 

Ft 


\p ( e+ -Uim + k B 




dp dpuj _ 
dt dxj 


(33) 

dpiii 

dt + 

d 

dxj 

[pUiUj + PSij - (r ij + T?/ s )] = 0 


(34) 

d 

+ dx j 

pu 3 

^ h + —UiUi + k' 9 '^j — Ui (jij + + ( q 3 + q B 9S ) 

= 0 

(35) 


where the sub-grid scale (SGS) heat flux, stress, and kinetic energy are defined as: 


q S j 9S = p(ujh-Ujti) 

T ij S = ~P (Fuj - UiUj) 

pk S9B = ^p{FFi ~ Fui) 


(36) 

(37) 

(38) 


These sub-grid scale terms must be modeled. However, unlike in the RANS equations where all of the 
turbulent effects were replaced by a turbulence model, the SGS terms only represent that portion of the 
turbulent motion that is at too small of a scale to be resolved by the computational scheme and mesh. 

A large-eddy simulation is a complex interaction between computational grid, numerical scheme, bound- 
ary conditions, sub-grid model, and filter. The key to the success of LES is properly resolving the large-scale 
turbulent structures that contain the majority of the turbulent energy in the flow. Maximizing the resolution 
of the simulation minimizes the effect of the sub-grid model. Because LES directly computes the large-scale 
structures, it better represents the turbulent flow physics and is not subject to the types of modeling errors 
associated with RANS simulations of jets and mixing layers. The effects of compressibility on turbulence 
and mixing will be largely captured directly, with only minimal errors occurring due to the sub-grid model. 
The primary LES issue specific to free shear flows is the proper resolution of the developing region. This re- 
quires sophisticated boundary conditions to simulate the initial turbulent inflow and high-resolution schemes 
together with high-quality grids to capture the small but important turbulent structures. 

A brief summary of the key LES issues is presented below. Georgiadis and DeBonis 103 include a more 
complete review of the subject. 


III. A. Numerical Scheme 

The numerical scheme used in LES/DNS calculations can have a significant effect on the computed results. 
In order to resolve the evolution of the turbulent structures without damping them, schemes that maximize 
wave resolution and minimize numerical dissipation are desirable. For this reason, high-order explicit finite- 
differences are often used to maximize the order of accuracy of the scheme. Dispersion relation preserving 
(DRP) schemes further optimize the discretization stencil coefficients to minimize the dispersion error over 
a limited range of scales. 110, 111 Compact schemes are also popular and increase the resolution for a given 
stencil size by implicitly computing the derivatives. 112 Very low dissipation methods based on kinetic energy 
conserving schemes have shown significant promise for low-order unstructured grid codes. 113 

Central difference formulations are desirable for these types of simulations because they lack inherent 
dissipation. However, they will be numerically unstable unless some form of artificial dissipation is added. 
This is frequently done through the application of a filter function which is designed to work with a particular 
scheme. Upwind schemes are naturally more dissipative than central difference methods and will damp the 
turbulent structures. Use of higher-order upwinding is usually necessary to allow the growth and evolution 
of the turbulence. Note that regardless of the scheme used, the order of accuracy only indicates how the 
error will be reduced with grid refinement and does not indicate the magnitude of the error for a given grid 
size. 

Compressibility introduces another complexity for large-eddy simulation, namely how to handle shock 
waves that are present in the flow. In central difference solvers, a trigger based on the pressure gradient is 
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often used to activate localized artificial dissipation or to switch to a blended central/upwind method in the 
region near a shock. Pure upwind-based methods tend to be more stable near shocks due to their higher 
levels of inherent dissipation. Total variation diminishing (TVD) limiters may also be used to reduce the 
pressure fluctuations on either side of the shock. Applying these fixes near shocks is a greater challenge in 
LES simulations than in RANS because any additional damping will act to weaken the resolved turbulent 
structures. Therefore, care must be taken to only apply enough dissipation to keep the scheme stable in 
the shock region. The localized artificial diffusivity (LAD) technique 114 provides shock stability without 
affecting turbulent structures by adding artificial dissipation through the bulk viscosity term, leaving the 
shear viscosity and the stress tensor unaffected. 

III.B. Computational Mesh 

Numerical schemes are sensitive to the computational mesh used. Each scheme has its own requirements 
for the number of points to resolve a particular wave number, grid skewness, and rate of grid stretching. 
Higher-order methods used in LES/DNS often have wide discretization stencils that introduce additional 
grid requirements and make the treatment of complex geometry challenging. Isotropic grids are best for 
capturing the rotating nature of the turbulent structures as they convect downstream. 

For LES/DNS, the range of scales that must be resolved is Reynolds number dependent. As the Reynolds 
number increases, the Kolmogorov length scale, which represents the smallest eddies, will decrease. While 
LES does not resolve all the way down to the Kolmogorov scale, the range of “important” scales will likewise 
grow. For higher Reynolds number flows, this means that a finer grid will be needed to resolve a similar 
percentage of the overall scales or else the sub-grid model will be relied upon to model their effect. 

In a time accurate simulation such as LES, the size of the step used to advance the solution in time 
is directly proportional to the smallest grid spacing in the domain through the Courant-Friedrichs-Lewy 
(CFL) stability condition for the numerical scheme used. Increasing the grid resolution not only increases 
the computational time per iteration to solve the larger grid, but it also increases the number time steps 
necessary to achieve the same total simulation time. It is important to note that the total simulation time 
needed is determined by the size and frequency of the large-scale structures, and remains fixed for a given 
flow. 

Wall bounded flows likewise present a challenge for LES, because of the range of scales involved. This 
also becomes important for free shear flow applications, because the state of the upstream boundary layers 
plays an important role in the development region of the shear layer. Even if the unsteady flow entering the 
free shear region is prescribed, the downstream grid must be fine enough to resolve the incoming turbulent 
eddies and give them a chance to grow. Otherwise the numerical dissipation will dampen if not destroy the 
turbulent fluctuations. 

As an example, consider the case of an axisymmetric jet exiting from a convergent nozzle from reference 
115. The flow conditions for this subsonic unheated jet are very similar to setpoint 7 in table 1. Figure 13, 
which compares the centerline velocity decay, reveals a potential core that is a little too short and a decay rate 
that is too fast in the region near the end of the core, but then a mixing rate that seems to improve farther 
downstream. Examination of the root mean square (RMS) velocities downstream of the nozzle lip (figure 14) 
reveal initial values that are much too high, but then provide better agreement farther downstream. Because 
the grid near the nozzle exit is not fine enough to capture the small-scale turbulence from the boundary 
layers, the LES simulation sees this as a more laminar-like jet. The elevated RMS velocities downstream of 
the nozzle lip (figure 15) are characteristic of laminar to turbulent transition. This laminar-like inflow results 
in the shorter potential core, but as more of the energy-containing eddies are resolved farther downstream 
the rate of centerline velocity decay improves. 

III.C. Boundary Conditions 

Boundary conditions present a particular challenge for LES/DNS. Rarely are the unsteady inflow conditions 
known to the level of fidelity needed, and computing the entire unsteady domain is not usually feasible. In- 
stead several techniques have been developed to provide approximate inflow conditions. Synthetic turbulence 
methods impose an unsteady velocity profile based on a specified mean flow and Reynolds stress profile. The 
specified profiles may come from experimental data or a RANS simulation that is run a priori or in a coupled 
hybrid RANS/LES fashion. Various techniques have been developed to compute the fluctuating quantities 
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Figure 13: LES axisymmetric jet centerline velocity from reference 115. 



Figure 14: LES axisymmetric jet lip line turbulence from reference 115. 
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(c) Streamwise RMS Velocity (LES) (d) Radial RMS Velocity (LES) 

Figure 15: Experimental and LES axisymmetric jet turbulence from reference 115. 


that reproduce correct Reynolds stresses: burst models, 116 synthetic eddy methods, 117 and digital filter- 
ing. 118 Another method, used with wall boundary layers, involves a recycling technique whereby a section of 
the upstream boundary layer is simulated and the outflow of that region is rescaled and re-used at the inflow 
boundary to create a realistic inflow condition to the domain of interest. 119, 120 Care must be taken to ensure 
that the upstream region is not so short as to artificially limit the size of turbulent structures it allows. The 
recycling method can also be difficult to implement when the flow is not predominantly one-dimensional or 
when complex geometries are involved. 

Reflecting waves from far-held and outflow boundaries can propagate back into the region of interest and 
contaminate a solution. Numerical boundary treatments constructed using characteristic theory can be used 
to reduce unwanted wave reflections. Alternatively, many simulations employ a “sponge layer” of grid cells to 
absorb the outgoing waves. By increasing the grid spacing in these outflow regions, the dissipation inherent 
in the numerical scheme will help to damp the outgoing waves. Sometimes additional explicit dissipation 
must be added along these boundaries. 

III.D. Sub-grid Modeling 

In LES, sub-grid scale models are used to account for the unresolved unsteady motion, which has a net 
dissipative effect on the resolved how. Explicit sub-grid models are formulated to dissipate the appropriate 
amount of energy from the large-scale turbulence based upon the resolved motion in the simulation and the 
filter width, which is typically related to the grid size. Meneveau 121 categorizes sub-grid models into three 
general categories: 1) eddy viscosity models, 2) Germano identity and dynamic models, and 3) similarity, 
nonlinear and deconvolution models. Eddy viscosity models, such as the original Smagorinsky 122 model or 
Vremans model, 123 are the most widely used models for free shear hows. Bogey and Bailly 124 showed that 
the additional viscosity imposed through the sub-grid model lowered the effective Reynolds number of a jet 
how, altering the turbulent structures and resultant jet noise. Sub-grid models based on categories 2 and 3, 
such as the dynamic Smagorinsky model 125 and the approximate deconvolution model 126 are typically less 
dissipative and have shown good success. 

A popular alternative to sub-grid modeling is implicit LES (ILES). ILES encompasses a wide variety of 
approaches that do not employ a traditional sub-grid model but instead rely upon the dissipation in the 
simulation’s numerics to dissipate the large-scale turbulence. A popular subset of ILES methods is monotone 
integrated large-eddy simulation (MILES) 127 which employs monotone numerical methods such as Flux 
Corrected Transport (FCT), Monotone Upstream-centered Schemes for Conservation Laws (MUSCL), and 
Piecewise Parabolic Methods (PPM). In these methods the energy dissipation is implicit in the numerical 
scheme. Other ILES methods use centered difference methods and must rely on an explicitly added artificial 
dissipation term, such as a low pass filter. 115 ' 128 Results with any implicit LES method are obviously very 
dependent on the numerical scheme, but there are many excellent examples in the literature. Comparisons 
of implicit and explicit modeling using the same numerical method and grid show that the implicit method 
is much less dissipative. 128 
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III.E. Filtering 


Filtering separates the flow into large-scale resolved structures, that are directly resolved by the code, and 
small-scale unresolved structures which are modeled. Applying a filter to the Navier-Stokes equations is 
fundamental to the derivation of the LES equations. However, the need to explicitly filter the solution as 
part of the simulation process is less clear. Many researchers use the implicit filtering approach, arguing 
that the discretization process acts as a filter and no further filtering is necessary. 83 Others explicitly apply 
a filter to the solution at each time step. Explicitly filtering the solution is the more rigorous approach and 
gives the user control over the filter function and width. But in many codes, the real benefit of filtering is 
to provide artificial dissipation for stability, since it removes the small unresolved (or under-resolved) waves 
that lead to instability. Successful simulations have been performed with both approaches. 

IV. Future Research and Validation Needs 

The turbulence modeling community has periodically come together to evaluate turbulence models on 
a common set of problems in order to assess the current state of the art. A few well known gatherings 
include the 1968 AFOSR-IFP-Stanford Conference on the Computation of Turbulent Boundary Layers (often 
referred to as Stanford Olympics I), the 1972 Conference on Free Turbulent Shear Flows, and the 1980-81 
AFOSR-HTTM-Stanford Conference on Complex Turbulent Flows (Stanford Olympics II). More recently, a 
number of AIAA workshops have been held on the topics of drag prediction, shock wave turbulent boundary 
interaction, and propulsion aerodynamics. These activities have helped to identify shortcomings of existing 
turbulence models as well as problems with the way certain turbulence models have been implemented in 
different codes, including in some cases additional model corrections that were used by default but not fully 
documented. There is now an AIAA Turbulence Model Benchmarking Working Group which has a website 
at http://turbmodels.larc.nasa.gov/ devoted to disclosing model formulations, implementation issues, 
and verification cases. Users and developers now have a centralized source for verifying that their code and 
turbulence model will produce standardized results. The authors own work with other developers and codes 
have identified implementation issues that are sometimes problem specific. Due to the complexity of modern 
CFD codes, this type of collaboration is critically important and should continue. 

One on-going problem with model development is locating quality data for use in calibration and valida- 
tion. At present, there is no centralized repository for experimental data. While a few publishers now support 
the archival of on-line data, the amount of information one can include in a journal article is quite limited. 
NASA and/or AIAA should consider how to bring experimentalists together to analyze data for a series of 
fundamental cases relevant to turbulence model development, evaluate data quality and completeness, and 
if possible arrive at a set of consensus data that the modeling community can easily access. 

Additional experimental data is needed to further improve understanding of turbulent mixing layers 
and jets and aid in model development. For incompressible mixing layers, additional investigations are 
needed to better understand the relationship between different mixing layer thickness definitions, the effect 
of density ratio on growth rate, and the variation of peak turbulence statistics in the self-similar region. 
For compressible mixing layers, additional data from multiple facilities are needed to assess the changes in 
turbulence anisotropy. Care should also be taken to gather the type of detailed data necessary for calibration 
and/or validation of future modeling techniques. This includes: controlled and well-documented inflow 
conditions including the turbulence state and statistics; examination of both the development and self- 
similar regions; and measurements beyond just the mean flow and shear stress. Future model development 
will require knowledge of the three components of velocity in order to evaluate the mean stress and strain 
rate tensors and each of the Reynolds stresses, as well as higher-order correlations that contribute to the 
Reynolds stress transport equation. 

LES and DNS offer the most promise for providing “complete” flowfield information for use in RANS 
model development. As computational resources continue to advance and these methods are applied more 
frequently and at higher Reynolds numbers, the simulations should also be used as a means for evaluating 
RANS model approximations. Single-point correlations could be computed for the higher-order moments 
that are not readily available experimentally. These correlations could then be used to examine the validity 
of various model approximations throughout the flowfield (i.e., the Boussinesq approximation that relates 
the stress and strain rate tensors, the weak equilibrium assumption used in the formulation of algebraic 
Reynolds stress models, modeled forms of the pressure-strain correlation tensor, various gradient diffusion 
approximations, etc.). 
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V. Conclusions 


Turbulent free shear flows, particularly mixing layers and jets, play a vital role in aeropropulsion appli- 
cations. Experimental observations, especially at low Reynolds number, indicate the presence of large-scale 
coherent structures whose formation is triggered by the growth of instabilities in the initial part of the shear 
layer. As the flow progresses downstream, secondary instabilities are triggered that lead to the formation 
of streamwise vortices. The interaction between these various turbulent structures cause a distortion of 
the flow, increased three-dimensionality, and vortex stretching that lead to the development of small scale 
turbulence. Far downstream, the flow reaches a fully developed state where the profiles of mean velocity and 
turbulence statistics achieve self-similarity when plotted using appropriately scaled variables. In this region, 
the shear layer grows linearly with downstream distance. 

Reynolds-averaged Navier-Stokes (RANS) methods are widely used in production CFD, and are fre- 
quently able to capture trends exhibited by experiments when the geometry or flow conditions are altered. 
However, the accuracy of existing turbulence models still leaves much to be desired. These models fail to 
capture a number of basic, yet critical, flow features including: differences between planar and axisymmetric 
shear flows, sensitivity to compressibility effects on the flow, and temperature effects on shear layer growth. 
Attempts to resolve some of these issues by adding explicit dilatation terms, vortex stretching terms, and 
temperature corrections have had only limited success. There is ample experimental evidence that many of 
these phenomena are accompanied by changes in the organization of the large-scale turbulent structures. As 
such, RANS modeling efforts will likely need to move beyond simple eddy viscosity models which primarily 
affect the flow through changes in the turbulent shear stress. Reynolds stress models might provide a better 
representation of these flows since they bypass the linear stress-strain relationship imposed by the Boussinesq 
approximation used in eddy viscosity models, allow for turbulence anisotropy, and account for changes in 
the pressure-strain correlation tensor. It might be possible to further extend these models by making the 
closure coefficients sensitive to local flow invariants. 

Large-eddy simulation (LES) methods do not rely upon the numerous modeling approximations used by 
RANS models. Instead, the unsteady large-scale turbulent motion is resolved on the computational mesh 
and only the smaller scale turbulence is modeled. These methods naturally retain more of the physics 
described by the Navier-Stokes equations. However, LES is still subject to some uncertainty in physical 
modeling. Correctly specifying the unsteady inflow conditions is crucial for free shear flow simulations. How 
those conditions should best be approximated from only mean or statistical data remains an unknown. The 
method used to model the sub-grid scale contributions is also very important, whether that is done implicitly 
by the inherent dissipation in the numerical scheme or through explicit model terms and coefficients. LES is 
also sensitive to a number of numerical issues including: the dissipation of the numerical scheme, the use of 
explicit filtering particularly on non-uniform grids, and the resolution and quality of the computational mesh. 
The high-order numerical schemes often used for LES do not lend themselves well to complex grid topologies 
or unstructured grids. Current limitations in computational resources restrict the practical viability of LES 
methods to low Reynolds number flows. For free shear flow applications, adequate resolution of the small- 
scale turbulence in the initial part of the shear layer remains a challenge. As advances in computational 
resources allow for higher resolution grids and higher Reynolds number simulations, LES should become a 
more widely used tool for CFD simulations. 
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